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Abstract 

We examine a first order differential equation with respect to time coming up in the description 
of magnetic islands in magnetically confined plasmas. The free parameters of this equation are 
obtained by employing Bayesian probability theory. Additionally a typical Bayesian change point 
is solved in the process of obtaining the data. 
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I. INTRODUCTION 



Magnetic islands are structures appearing on resonant surfaces of plasmas in toroidal 
magnetic confinement devices. They have been found to limit the maximum achievable 
energy which can be stored in a fusion plasma and may therefore be a problem for a future 
reactor. Concepts of stabilizing the plasma in order to handle these instabilities include 
electron cyclotron current drive which can only be useful if it is accurately adjusted to the 
needed quantity. Therefore a thorough understanding of the island is necessary. 

The time dependence of the magnetic island width W is theoretically described by the 
generalized Rutherford equation [Q]. This first order nonlinear differential equation with 
respect to time contains in our case three free parameters which have to be determined from 
measured data since theoretical considerations can only provide estimates for these values. 
They are assigned to three terms describing stabilizing and destabilizing effects in the plasma, 
i.e. the bootstrap effect (with parameter Qbs), the Glasser-Greene- Johnson effect {accj) and 
the polarization currents (apoi). We use a simplified form of the Rutherford equation which 
comprises the relevant dependencies on the parameters = (a^s, accj, dpoi) only. The full 
account of all physical constants and time dependent quantities may be found in [Q. 

dW(t) W(t) ,,1 ,,1 

= const + «^^^^^(^) ll^-p^ - ^GGjCGGjit)^^ - apolCpolit)^, . (1) 

W^mm=l-8cm is the minimum width of an island. The variables cbs, cggj and Cpoi contain 
fundamental constants and time dependent input quantities like plasma temperature or 
pressure. 



II. GENERATING THE DATA 



The data is obtained from the so called Mirnov coils which are distributed poloidally 
around the torus and measure any change of the poloidal magnetic field. The time variation 
dm/dt of the magnetic flux m{t) through the Mirnov coil is proportional to the recorded 
signal. We are interested in the time evolution of the amplitude m of the integrated signal 
which can be connected to the magnetic island width W via 

Wit) = y!!!W^ , (2) 
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FIG. 1: Signal from Mirnov coils for discharge ^^12238 of the plasma device ASDEX Upgrade for 
different times. At ~ 2.6s the island has not formed yet, while at ~ 4.6s one can see the induction 
of the magnetic island in the coil signal. Finally at ~ 5.6s the island has disappeared. The peaky 
structures have to be removed from the signal for further process. 

where rrio is the offset of the magnetic signal and b a proportionality constant. Additional 
information about the absolute size of the magnetic island for a certain time comes from the 
electron cyclotron emission (ECE) diagnostic. From this we know that Wece = 7cm within 
a range of /S.ece = 1cm. This information will be used later in setting up a prior. 



A. Extracting the data from the Mirnov signal 

The original signal from the Mirnov coils is shown in Fig. |T[ A closer look (upper graph in 
Fig. 1^) reveals two kinds of structures: On the one hand peaks at intervals of approximately 
3-5ms which are due to an edge plasma phenomena where energy and particles are expelled 
out of the confined region, and on the other hand the signal originating from the change 
of the magnetic field which shows sinusoidal behavior (interval approximately 0.08ms) with 
an amplitude connected to the magnetic island width. This is the information we want to 
extract. First one has to identify the positions of the peaky structures. Since the height 
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FIG. 2: Original signal (upper graph) and after Fourier transformation back and forth where the 
higher frequency part was removed by filtering (lower graph). The peaky structures are easily 
identified (full circles). 

of the peaks is not everywhere larger than the highest amplitude of the sinusoidal signal 
we can not simply look for all points which are higher than a certain level. Fortunately 
the two structures live on two different scales in frequency domain. Therefore we Fourier 
transform (FFT) the complete data set and discard all the higher frequencies which refer to 
the sinusoidal structure (see Fig. Q). Back transformation gives then a signal where peaks 
are easy to identify. With the peak positions at hand we are set to go for the amplitude 
of the sinusoidal structure in between two peaks. Again Fourier transformation is applied 
where in addition we integrate over time and are finally left with the magnetic signal m 
shown in Figs. |] and § . 

B. Finding the valid range of the model 

The Rutherford equation (|I|) describes the dynamics of a magnetic island considering 
certain plasma physics effects. However, at the onset of the mode the magnetic island 
is not stabilized and subjected to fluctuations which are not covered by the model used. 
We therefore have to identify the region in which the Rutherford equation is valid. Fig. |^ 
depicts the amplitude of the magnetic signal. The left part differs from the right one where 
the island has stabilized in amplitude and noise and we have to look for the certain time 
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FIG. 3: Amplitude of the magnetic signal of discharge #12238. After the onset of the magnetic 
island the signal has not saturated and fluctuates (left side). The analysis of the Rutherford 
equation can only be performed for the stabilized part on the right side. The full line is the 
posterior distribution for the change point Nc- The dashed line is a fourth order polynomial in the 
respective regions. 



incident tc = t{Nc) when the change of the behavior happens - a typical Bayesian change 
point problem. Since we are out for the change point only we describe the time variation of 
the amplitude of the magnetic signal by low order polynomials 
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For convenience we use matrix notation in the following, with matrix elements {M}jt_i j = 
and cF = {ao, ax)- The index < (>) denotes time points before (after) the change point. 
The data is corrupted by noise: 



m> = M>/3 + £> . 



(5) 
(6) 



Again we assume (£</>) = and (£:</>)^ = cr|y^. Then the hkehhood reads 

1 1 
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•exp |-^(m< - M<«)^(m< - M<a) - ^(m> - M>/3)^(m> - M>/3)| . (7) 

We need the posterior distribution for the change point. With the help of Bayes theorem 
we get 

pmm,I) = ^^p{m\N,,I) . (8) 

The nominator in the fraction is the prior distribution in absence of any data which is a 
constant p{Nc\I) = const since no change point is preferred, but hmits the possible values to 
K+l<Nc<N — K — 1. The marginal likelihood p{Tn\Nc, I) is obtained by marginalizing 
over all parameters in (|^. 

p{m\Nc,I) = j da d(3 da<: da> p{m\Nc,d,f3,a<:,a>,I)p{d,/3\I)p{a^,a:^\I) . (9) 

For the prior in a and (3 we take a constant but use Jeffreys prior for p(cr<, cr> |/) = 1/ (cr<cr>). 
All integrations can be performed analytically and finally yield 

1 r(M) 



p{Nc\rn, I) oc 



detM^M< [m^m<-m5M<(M^M<)-iM^m<] 

^ 2 ^ (10) 



MetM^M> [m^m>-m^M>(M^M>)-iM^m>] 

The posterior change point distribution is shown in Fig. ^for a polynomial with order K = 4. 
We checked order three to five to get the same result. 

III. LIKELIHOOD AND PRIOR 

The measured quantity in Eq. (0) is the magnetic signal m with measurement uncertainty 
e. It is given in the form of a time series with successive events, where we can write 

mi = mo + b-Wi{ay + e , i = l,...,N . (11) 

Assuming that (e) = and (e^) = cr^ we get by virtue of the principle of maximum entropy 
the likelihood § 

p{m\mo, b, a, a, I) = (^2t^^2^n/2 Y1 K ~ - &W^(a)y | . (12) 



Next step is the assignment of prior distributions for the conditional dependencies on rrio, b, 
the free parameters a and a. Due to the above mentioned treatment of Fourier transforming 
the data back and forth we loose any information about the actual scatter originating from 
the measurement process. All we know is that a variance a exists and that it functions like 
a scale parameter - justified reasons for employing Jeffreys prior. 

p{a\I) = - . (13) 
cr 

Looking at the data before or after the formation of the magnetic island provides an estimate 
rho and its uncertainty for the offset and leads to a Gaussian prior distribution 

p{mo\mo,arh,,I) = exp <j ^ '> . (14) 



The ECE measurement is a constraint on the range of the proportionality constant h. Given 
a certain value for the offset rrio we insert Wece in Eq- (0) which gives an estimate bECE = 
{riiECE — ^o)/Wece- The ECE measurement uncertainty Aece may be used in order to 
set up an upper and lower limit: b^p/iow = {^ece ~ i^o) / {Wece — / + ^ece)"^ A convenient 
prior function given an estimate within boundaries is the beta prior . However, it operates 
for values between and 1 only, so we have to renormalize x{b) = {b — biow)/ibup — how)'- 

p{b\mo,WECE,AECE,mECE,I) = YTT^ " ^^""'[1 " <b)Y~' , (15) 

where 

(1 - /ife)/ib - 0-^2 l-/ib 
U = 2 ■ fl'b , V = U (16) 

and 

Hb = x{bECE) , cTfe = - — ^- — Aece ■ (17) 

Vup — Olow yy ECE 

From theoretical considerations we have some idea about the quantities of the free parame- 
ters but unfortunately only for certain ideal configurations of the confined plasma, ass = 1-7, 
o-GGJ = 6-5/9 and apoi = 7 are provided by literature ^ . The maximum entropy principle 
gives us in this case an exponential function. 

p(a|ao,J) = J] — exp|-^| . (18) 
j=i ^ J 



IV. PARAMETER ESTIMATION 



We are out for the parameters of the Rutherford equation a together with ttIq and h 
from Eq. (j^). In the expectation value for a component oi 9^ = {rrio^b, Ubs, clggj, o,poi) we 
marginahze over all variables entering the likelihood Eq. (^2|): 



SA0jp(e,\in,I) fdS Jdap(e,a\m,I) J ' 

p{9) may be used as a sampling density in Markov chain Monte Carlo (MCMC). Invoking 
Bayes theorem 

p[e, a m, I) = -— — 20) 

p[m\Ij 



gives 
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The full prior in Eq. (^) disentangles into the functions given in Eqn. (|l^ , [l5|Jl8| , p^) 

p{0,a\I) = p{mo\fno,afhoJ)p{^'m.o,WECE, ^ECE,mECEj)p{a\aoJ)p{^^^ . (22) 
The integration over a can be treated analytically and results in 
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/ do" p{m\6, a, I)p{a\I) oc < [mj — rrio — 6VFi(a)^] 

I i=i 



(23) 



The final integrations over the parameters 9 are performed numerically by employing 
MCMC, while the first order differential equation (P is solved applying second order Runge- 
Kutta method. 



V. RESULTS 

The analysis is performed for discharge #12238 of the plasma device ASDEX Upgrade. 
Fig. § depicts the dynamics of the magnetic island width. The thin line is the width 
obtained from the magnetic signal m employing Eq. ^ with the expectation values of nio 
and h from the analysis (see table |). Only that time interval of the complete signal is 
examined which comprises the island after is has stabilized until the temperature signal 
shows decoupling from the behavior of the collapsing island. The comparison with the 
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FIG. 4: Magnetic island width for discharge #12238 of ASDEX Upgrade. 

experimental data (thin line) gives a very good agreement. The accompanying parameters 
are given in table |. 



VI. SUMMARY 

Bayesian analysis was employed in order to identify the valid region in a data set needed 
for further examinations. The evolving data was used to determine free parameters in the 
Rutherford equation, a first order nonlinear differential equation describing the magnetic 

mo [10"^] b ass acGJ OpoZ 

0.582 ± 0.075 0.184 ± 0.004 0.770 ± 0.007 1.26 ± 0.04 0.783 ± 0.035 

TABLE I: Expectation values with error margins for the magnetic offset mo, Unear factor h and 
the three parameters of the Rutherford equation. 



island dynamics in toroidally confined plasmas. 
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